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Abstract 

The geological record documents links between glacial cycles and volcanic productivity, 
both subaerially and, tentatively, at mid-ocean ridges. Sea-level-driven pressure changes 
could also affect chemical properties of mid-ocean ridge volcanism. We consider how 
changing sea-level could alter the CO 2 emissions rate from mid-ocean ridges on both the 
segment and global scale. We develop a simplified transport model for a highly incompat¬ 
ible trace element moving through a homogenous mantle; variations in the concentration 
and the emission rate of the element are the result of changes in the depth of first silicate 
melting. The model predicts an average global mid-ocean ridge CO 2 emissions rate of 
53 Mt/yr or 91 Mt/yr for an average source mantle CO 2 concentration of 125 or 215 ppm 
by weight, in line with other estimates. We show that falling sea level would cause an 
increase in ridge CO 2 emissions about 100 kyrs after the causative sea level change. The 
lag and amplitude of the response are sensitive to mantle permeability and plate spread¬ 
ing rate. For a reconstructed sea-level time series of the past million years, we predict 
variations of up to 12% in global mid-ocean ridge CO 2 emissions. 


1 Introduction 

Glacial cycles transfer ~ 5 x 10 19 kg of water between the oceans and ice sheets [Tushingham 
and Peltier, 1991], leading to accumulation and ablation of kilometres of ice on the continents 
and sea-level change of ~ 100 m. I 11 Iceland, for example, the pressure change associated with 
melting of the ice sheet since the last glacial maximum had well-documented consequences for 
the volcanic activity [Sigvaldason et ah, 1992, Jull and McKenzie, 1996, Maclennan et ah, 2002] 
and lava geochemistry [Maclennan et ah, 2002]. More broadly, continental volcanism in both the 
northern and southern hemispheres shows increased activity associated with the last deglacia¬ 
tion [Gardeweg et ah, 1998, Jcllinek et ah, 2004, Huybers and Langmuir, 2009]. Huybers and 
Langmuir [2009] and Lund and Asimow [2011] hypothesised that the pressure variations caused 
by changing sea level during glacial cycles would affect mid-ocean ridge (MOR) volcanism. 
Crowley et ah [2015] documented variations of bathymetry near the Australian-Antarctic ridge 
that are possible evidence of such glacial effects. 

A simple argument shows that variations in crustal thickness and sea-floor relief should be 
expected to result from sea-level variation. The melting rate of a parcel of mantle beneath a 
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MOR is proportional to its depressurisation rate. As the parcel upwells, it depressurises due 
to the decreasing height of rock above it. The rate of change of pressure due to upwelling is 
the gravitational acceleration g times the mantle density times the upwelling rate (~ 3 cm/yr). 
Sea-level variation can modify this depressurisation rate: the pressure change due to varying 
sea level is the product of g, water density and the rate of change of sea level (up to 100 m in 
10 kyr or 1 cm/yr). Water is about one third the density of the mantle and sea-level changes 
can be equivalent to one third the mantle upwelling rate, implying that sea-level changes can 
modify depressurisation rates, and hence melting rates, by up to ±10%. Crowley et al. [2015] 
apply a paleo-sea-level reconstruction to a simulation of MOR melting and melt transport, 
with melting rates varying according to the variable pressure exerted by sea level. This leads 
to varying melt flux at the ridge, predicting variations in crustal thickness consistent with 
bathymetric observations of sea-floor relief. 

Given this evidence for glacial cycles affecting MOR melt production, it is reasonable to 
consider if, as in Iceland [Maclennan et ah, 2002], the chemistry of the lavas is also affected. To 
investigate this we develop a model of the transport of a highly incompatible element from the 
asthenosphere through the melting region to the MOR. Highly incompatible elements partition 
strongly into the melt, rather than remaining in the residual solid. Approximating this as 
complete incompatibility creates useful simplifications in modelling. For instance, a completely 
incompatible element’s path through the melting region is entirely determined by the motion of 
the melt, without any need to consider that element in the solid or partitioning between phases. 
Furthermore, for small perturbations to the melting rate (such as those caused by sea-level 
change) the mass-flow rate of the element through the MOR is constant (see appendix A.4). 

In a simple model of melting beneath a mid-ocean ridge, a parcel of mantle upwells adia- 
batically beneath the ridge axis, cooling slightly due to its expansion. The pressure-dependent 
solidus temperature of that parcel decreases as it ascends; at some depth (or, equivalently, pres¬ 
sure), the temperature of the parcel is equal to its solidus temperature. This depth, thought 
to be around 60 km, is called the depth of first silicate melting. With further upwelling, the 
parcel’s temperature exceeds the solidus and it partially melts. As soon as the first increment 
of melt is produced, 100% of the completely incompatible element that was locally present in 
the solid mantle is transferred to the melt. Because the mantle is permeable and the melt is less 
dense than the residue, the melt ascends faster than the solid, segregating from its source. Melt 
segregation and transport of incompatible elements thus begins at a pressure-dependent depth. 
More specifically, melt segregation begins at a fixed pressure, but the depth corresponding to 
this pressure can change. We assume that there is no isostatic rebound associated with sea-level 
change and hence mantle upwelling is constant. Variations in sea level will therefore cause the 
depth of first silicate melting (and initiation of melt segregation) to rise and fall. The rate at 
which mantle crosses this boundary and delivers its content of incompatible elements to the 
melting region is the mantle upwelling rate minus the rate of upward motion of the boundary. 
So in this model, variations in the depth of first silicate melting cause variations in the flux 
of an incompatible element. As sea level falls, the depth of first melting increases, upwelling 
mantle crosses into the melting region faster, and the flux of incompatible element increases; 
the reverse is true for sea-level rise. Any perturbation to the melting rate within the melting 
region does not alter the mass of the incompatible element in the melt, it only dilutes (or con¬ 
centrates) the element. Variations in melt-transport rate associated with melting perturbations 
are a secondary effect and are not considered in detail here (though see appendix A.4). 

A more nuanced view of melting may disagree with this simple story in some of the details, 
especially with the inclusion of volatile elements that are present in small concentrations in the 
mantle. Experimental evidence suggests that CCH-rich melt forms at ~ 250 km depth [Dasgupta 
et ah, 2013] and has a low viscosity that rises sharply with silica content at shallower depth 
[Kono et ah, 2014], If such carbon-rich melts can segregate from the solid mantle it would 
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complicate the role of the transition to silicate melting at around 60 km. However, it remains 
an open question whether oxygen fugacity allows such melts to form and, if they do form, 
whether such tiny melt fractions can segregate from the solid mantle. Dasgupta et al. [2013] 
suggests carbonatite melt fractions reach ~ 0.03 wt% deep in the mantle below ridges, which 
is at the lowest limit of carbonatite melt interconnectivity of 0.03-0.07 wt% in ~ 0.05 mm 
olivine grains [Minarik and Watson, 1995]. The additional presence of water might increase 
the melt fraction to 0.06-0.1 wt% [Dasgupta et ah, 2013] by 150 km depth, but the threshold 
for interconnectivity of such melts is not known. In our calculations, we assume that these 
melt fractions do not segregate from the solid mantle until the onset of volatile-free peridotite 
melting at ~ 60 km raises the melt fraction, creating an interconnected, permeable network of 
pores. 

Among the highly incompatible elements, we focus on carbon despite its active role in the 
thermodynamics of melt production because variations in CO 2 emissions from the solid Earth 
are potentially important to understanding past variation of the climate. The solid Earth con¬ 
tains 10 1 °-10 11 Mt carbon [Dasgupta and Hirschmann, 2010], orders of magnitude more than the 
atmosphere (0.6 x 10 6 MtC [Solomon et ah, 2007]) and the oceans (4 x 10 7 MtC [Solomon et ah, 
2007]). Solid-Earth carbon emissions from MORs are estimated as ~ 25 MtC/yr [Coltice et ah, 
2004, Cartigny et ah, 2008, Marty and Tolstikhin, 1998] and from arc volcanoes ~ 20 MtC/yr 
[Coltice et ah, 2004], As the largest reservoir, the solid Earth’s carbon budget is known to con¬ 
trol atmospheric carbon on multi-Myr timescales; geological ages show a correlation between 
volcanic activity and atmospheric CO 2 concentration [Budyko et ah, 1987]. Hence there is 
evidence for both MOR volcanism being affected by glacial cycles and for the effect of volcanic 
CO 2 emissions on atmospheric CO 2 concentration. While we focus on CO 2 in our model, the 
same theory applies equally to other highly incompatible elements such as U, Th, Nb, Ba, and 
Rb. 

The model is developed under the guiding principle that it should be simple enough that 
the connections between the assumptions and the outputs are readily traceable. The full model 
is assembled from independent, decoupled parts that capture the key physics with minimal 
complexity. Mantle flow is modelled by the passively-driven corner-flow solution [Batchelor, 
1967, Spiegelman and McKenzie, 1987]; lithospheric temperature structure and thickness is 
computed with a half-space cooling model [Turcotte and Schubert, 2002], To quantify melt 
generation and transport we use one-dimensional compaction columns [Ribe, 1985, Hewitt, 
2010] that are based on conservation of energy, mass, and momentum at steady state in a 
homogenous mantle. The outline of the melting region is given by a parameterised solidus 
[Katz, 2008]. A focusing width is applied such that melt focused to the ridge produces a 
maximum crustal thickness of 7 km. A detailed discussion of the assumptions made in deriving 
the model is presented below. 

To summarise the results, the model predicts that a section of MOR spreading at 3 cm/yr 
half-rate will see a change in the rate of efflux of highly incompatible elements ( e.g ., CO 2 ) of 
~ 8% for a linear sea-level change of 100 m in 10 kyrs. For reconstructed sea level data and 
the present distribution of plate spreading rates, the model predicts global MOR emissions to 
deviate from the mean by up to ±6%. These results are sensitive to the permeability of the 
mantle, which is a primary control on the rate of melt transport. There are good constraints on 
how permeability scales with porosity, but its absolute value at a reference porosity (1% here) 
is uncertain. We consider a broad range of values that includes the most extreme estimates. 

Section 2 details the model used to predict C0 2 emissions for a section of mid-ocean ridge; 
parameter values are stated in table 1. The behaviour of the model is demonstrated for simple 
scenarios of sea-level variation in sections 3.1-3.3; the model is applied to the global MOR 
system under a reconstructed sea-level history in sections 3.4 and 3.5. The results are discussed 
in section 4 and the key conclusions stated in section 5. 
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2 The Model 


Our goal is to develop a method to compute the C0 2 emission rate E C q 2 (mass per unit time 
per unit length of ridge) from a segment of ridge. To achieve this we require a model of C0 2 flux 
into the melting region and also of its transport to the ridge. We approximate the behaviour of 
C0 2 as perfectly incompatible, and hence there is no exchange of C0 2 between phases during 
melt transport. The rate of ridge emission of C0 2 can then be quantified by integrating the 
mass flux into the base of the melting region fco‘2 (mass/area/time) and using the travel time 
from the base of the melting region to the ridge. This is formulated as 

E co 2 i.t) = ‘ 2 [ fco 2 {t s ,x,U 0 )dx , (1) 

Jo 

where x is the horizontal distance from the ridge axis, Xj is the maximum distance over which 
melt is focused to the ridge axis, Uq is the half-spreading rate, and t is time. A parcel of 
melt arriving at the ridge axis at time t was produced by mantle that crossed into the melting 
region at time t 3 . The travel time of melt from the base of the melting region to the ridge is 
represented as r, which varies with lateral distance x from the ridge axis. Hence the source 
time is t s (x) = t — t(x). The factor of two in eqn. (1) arises from the symmetry of the melting 
region across the ridge axis. A sketch of half the melting region is shown in figure 1. 


x = 0 

z = 0 - 1 - 



W m (x) 


Figure 1: Sketch of the melting region. Two example melt streamlines are shown in green. 
Mantle upwelling rate into the base of the melting region W m is represented by grey arrows, 
the size of which shows the decreasing magnitude of W m with distance from the ridge axis. 
The red arrow indicates the variability of z m with respect to time, expressed in equation (6) 
(all other boundaries are steady state). The melt streamlines on either side of the maximum 
focusing distance Xf show melt flow to the ridge or frozen into the lithosphere. 


Equation (1) requires an expression for fco2- This is a product of the rate at which mantle 
material crosses the depth of first silicate melting z m and the C0 2 concentration in that mate¬ 
rial. For generality, we consider an expression that allows for volatile-enriched insipient partial 
melting beneath the depth of first silicate melting, though we will later exclude this scenario 
from consideration. Thus, fco 2 is written as 

fco,(t„x,Uo) = (w m {x,Uo) - (1-0)CS£+ (w m - i ^^)<pcggl , ( 2 ) 
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where W m (x, Uq) is the upwelling rate of the mantle, w m is the upwelling rate of insipient melt, 
Cqo 2 i s a mass concentration of CO 2 , and </> is the volume fraction of melt; all of these are 
evaluated at depth z m and distance x from the ridge axis. The first term in parentheses on 
the right-hand side of the equation is the rate at which mantle crosses into the melting region. 
The concentrations and melt fraction may be considered steady-state, constant values as long 
as w m > W m > ma x(z m (t)). This ensures that the rate at which material crosses the depth 
of first silicate melting is always positive or zero. Even the fastest sea-level changes on record, 
meltwater pulses during the last deglaciation, satisfy these conditions for most MORs: 20 m 
sea level change in 500 years gives z rn = 1.3 cm/yr. However, meltwater pulse events are not 
resolved in the reconstructed sea-level series that we consider in this paper [Siddall et al., 2010], 
so the conditions w m > W m > max(i m (f)) are true with only occasional exceptions for the 
slowest spreading ridges. With these conditions satisfied and assuming that either 4>(z m ) = 0 
or that w m = W m , equation (2) can be simplified to 

fco 2 (t s ,x,U 0 ) = (w m (x,U 0 ) - C c02 . (3) 

Here we should interpret Cqo 2 as the ma ss concentration of C0 2 in the solid mantle plus co¬ 
moving insipient melt, if it is present. This model could be modified to accommodate more 
complicated situations, but this is not given further consideration at present. 

The flux of C0 2 in equation (3) depends on the solid mantle upwelling rate at the depth of 
first silicate melting. Approximating this as passive (plate-driven) flow of isoviscous rock, man¬ 
tle upwelling is given by the corner flow solution [Batchelor, 1967, Spiegelman and McKenzie, 
1987]. The vertical component of this solution, evaluated at z = z m , can be written as 


W m (x,U 0 ) 


2U, 


0 


7 r — 2 a r — sin 2a r 




(4) 


where the lithosphere is represented as a wedge with angle a c to the horizontal. We follow 
Spiegelman and McKenzie [1987] in computing the wedge angle to approximately match the 
plate thickness at a specified distance from the axis. We use a c = tsxT 1 {zi/x w ), such that the 
wedge intersects the upper boundary of the melting region zi at the maximum width of melting 
region x w (see appendix A.3 for details). This definition of the wedge angle ensures that the 
lithospheric wedge does not overlap with the melting region, and that the upwelling rate is 
small but non-zero at the extreme width of the melting region (0 < W rn (x w ) <C Uq). 

An expression for the depth of first melting z rn is needed in equations (3) and (4). Taking 
decompression as the only influence on local mantle temperature prior to melting, we model 
the depth of first melting as the intersection between an adiabatic temperature profile and 
the solidus temperature profile. We approximate both profiles as linear with respect to depth 
(details in appendix A.l) to obtain 


z 


m 


T Ts ie f A ^ Pw_g 

ipg-^f) p 


(5) 


where T is the mantle potential temperature, Ts ref is the solidus temperature at reference mantle 
composition and surface pressure, 7 is the Clausius-Clapeyron slope for the mantle, p is the 
mantle density, a is the coefficient of thermal expansion, c is the specific heat capacity, S is the 
sea-level deviation from a long-term mean, and p w is the density of water. The first term in 
equation (5) is the dry peridotite melting depth of ~ 60 km, and the second term is the shift 
in melting depth due to sea-level. The only time dependent variable in equation (5) is the sea 
level S, so differentiating gives 

d Z m Pw dS* , . 

= jm' (6) 
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These equations state that silicate melting begins at a fixed pressure, but the depth correspond¬ 
ing to this pressure varies as sea level rises and falls. 

The source time t s used in equations (1) and (3) depends on the travel time of melt to the 
ridge from any point x along the base of the melting region. For simplicity, we consider melt flow 
as following a vertical path from the base to the top of the melting region, then following a high 
porosity channel along the impermeable top of the melting region to the ridge axis, as illustrated 
by streamlines T\ and T 2 in figure 1. This is a reasonable approximation of numerically modelled 
streamlines for homogenous mantle [ e.g ., Katz, 2008], where buoyancy forces drive vertical fluid 
flow in the majority of the melting region, with the compaction pressure only becoming large 
enough to deflect melt flow from the vertical within a few kilometres of melt-impermeable 
boundaries [Sparks and Parmentier, 1991]. We consider flow along the high porosity channel 
as instantaneous, motivated by the high flow rates expected there, as compared to vertical flow 
rates in the rest of the melting region [Katz, 2008]. To compute the travel time from the base of 
the melting region to the base of the lithosphere, we use a ID compaction column from Hewitt 
[2010]. This model assumes Darcy flow and thermodynamic equilibrium for a two-component, 
homogenous mantle with a constant Clausius-Clapeyron slope. Following the reduced model of 
Crowley et al. [2015], we assume that small variation in melting rates due to sea-level change 
do not significantly affect melt velocity or travel time. Furthermore, in computing r we take 
z m as constant, because changes in z m due to sea level are only tens of metres, changing r by 
<C 1%. This gives a travel time 


7" (IKm, Zlj Z m ) 


( Vf V ( p 

\K 0 Apg) VniF m 



(7) 


where iy is the mantle melt viscosity, K 0 is the permeability of the mantle at the 1% reference 
porosity d 0 , Ap is the density difference between the solid and melt, n is the porosity exponent 
in the permeability relation (n ~ 2.7 ~ 3, Miller et al. [2014]), zi(x,Uo ) is the depth of the 
upper boundary of the melting region, and n is the adiabatic melt productivity (kg of melt 
produced per m 3 per m upwelling). This productivity is given by Hewitt [2010] as a ratio of 
thermodynamic parameters, and we retain the same parameter values here. 

To close the model we need an expression for the upper boundary of the melting region 
Zi(x). This boundary is located where the local temperature equals the solidus temperature, 
which is controlled by the thickness of the conductively cooled boundary layer that forms the 
lithosphere. Crowley et al. [2015] also included the effects of adiabatic decompression and latent 
heat removal, but this leads to an expression for xi(z) that cannot be inverted. However, as 
shown in appendix A.2, the depth of this melting boundary is approximated by an isotherm 
of the half-space cooling model, which can be expressed as Zi(x). This isotherm is hotter than 
the low-pressure solidus, but the temperature difference compensates for the change in solidus 
temperature with respect to depth. We use 


Zlix ' u ° ) = 2 \lw 0 erfc - 1 (|r^) ■ (8) 

where T) is the temperature of the upper boundary of melting region, assumed to be constant. 
To is the temperature of the ocean floor, and k is the thermal diffusivity of the mantle. Note that 
although we have neglected adiabatic decompression and latent heat of melting in equation (8), 
these are accounted for in the melting calculations. 


3 Results 

Below we demonstrate the behaviour of the model by a series of examples, then consider global 
CO 2 emissions for reconstructed sea level. We begin by modelling a unit length of mid-ocean 
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Parameter 

Value 


Parameter description 

c 

1200 

J/kg/K 

Specific heat 

Cco 2 

0.375 

kg/m 3 

CO 2 mass per m 3 of mantle (derived from 125ppm by weight) 

9 

10 

m/s 2 

Gravitational acceleration 

K 0 

10 - 13 - 10 " n 

m 2 

Permeability at 1% porosity 

K 

Ao( 0 / 0 o) n 

nr 

Permeability 

L 

4 x 10 5 

J/kg 

Latent heat of the mantle 

I'MOR 

61000 

km 

Total length of the mid-ocean ridge system 

n 

3 

- 

Porosity exponent in the permeability relation 

S 

[various] 

m 

Sea level, expressed as deviation from long term mean 

S 

[various] 

cm / yr 

Rate of change of sea level with respect to time 

T 

1648 

K 

Potential temperature of upwelling mantle 

To 

273 

K 

Ocean floor temperature 

Ts ie{ 

1554 

K 

Reference solidus temperature 

U 0 

< 8 

cm / yr 

Plate half-spreading rate 

x f 

10-70 

km 

Width of region, at z m , from which melt is focused to the ridge 

x w 

10-350 

km 

Width of melting region 

Zm. 

60 

km 

Depth of first melting 

a 

3 x 1(T 5 

K - 1 

Thermal expansion coefficient for the mantle 

7 

60 x 10“ 9 

K/Pa 

Clausius-Clapeyron slope of the mantle 

Vf 

1 

Pas 

Mantle melt viscosity 

K 

10“ 6 

m 2 /s 

Thermal diffusivity 

n 

0.01 

kg/m 4 

Adiabatic melt productivity, kg of melt per m 3 of mantle 
per metre of upwelling 

p 

3300 

kg/m 3 

Mantle density 

Pc 

2900 

kg/m 3 

Oceanic crust (mean) density 

Pw 

1000 

kg/m 3 

Freshwater density 

Ap 

500 

kg/m 3 

Density difference between liquid and solid mantle 

00 

0.01 

- 

Reference porosity (volume fraction) 


Table 1: Parameter values for calculations. 


ridge in the absence of sea-level change. This defines the baseline state of the model that 
we compare against when applying sea-level forcing. The first dynamic example is a single, 
linear sea-level change, which illustrates key characteristics of the model and emphasises the 
importance of the melt travel time. The emissions curve for this example approximates the 
Green’s function for the model: the response to a step-change in sea level. We next consider 
simple, periodic sea level curves, demonstrating the system’s response to an oscillatory forcing 
as a function of the frequency of that forcing. We then compute the model’s response to a 
reconstructed Pleistocene sea-level record. 

Building on these calculations for a single section of mid-ocean ridge, we calculate predictions 
for a global model, created by summing sections with appropriate spreading rates over the 
full MOR system. We demonstrate this composite model by considering the global emissions 
response to a single sea-level change and then to the reconstructed Pleistocene sea-level curve. 

3.1 Constant Sea Level &; Baseline Emissions 

Baseline emissions are defined as the steady-state emission rate of CO 2 from the MOR for con¬ 
stant sea level. This represents the background state, which is disturbed after sea-level change. 
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Figure 2 shows baseline emissions as a function of plate spreading rate and demonstrates that, 
in the absence of changing sea level, CO 2 emissions per metre of ridge are approximately pro¬ 
portional to the half-spreading rate Uq. A faster spreading rate drives faster mantle upwelling, 
bringing more CO 2 into the melting region. Faster spreading rates also increase the width of 
the melting region (see appendix A. 3), leading to more melt and more CO 2 being focused to 
the ridge axis. However, not all melt produced is focused to the MOR; at some lateral distance, 
melts are frozen back into the lithosphere, rather than travelling to the ridge axis [e.g., Katz, 
2008]. We incorporate this detail by enforcing a maximum focusing width Xf such that melt fo¬ 
cused to the MOR will produce crust < 7 km thick in steady state, consistent with observations 
[White et al., 1992, Bown and White, 1994], Crustal thickness is calculated as the volume flow 
rate of melt to the ridge axis (kg m~ 3 per metre along MOR axis) times the density ratio of 
basaltic melt to oceanic crust, divided by the half-spreading rate. At half-rates Uq > 1 cm/yr, 
the focusing width is smaller than the width of the melting region. With the imposition of this 
limit on melt focusing, there is a slight change in the slope of the baseline emission curve at 
Uq = 1 cm/yr (fig. 2). 

For the range of Uq on Earth of up to 8 cm/yr, figure 2(a) shows baseline emissions of 
up to 2300 kg m _1 yr -1 . This is calculated using 0.38 kg CO 2 per m 3 of upwelling mantle 
(125 ppm CO 2 by weight, hereafter ppmw), based on an average of MOR source mantle of 
50-200 ppmw [Dasgupta, 2013, Cartigny et al., 2008, Saal et ah, 2002, Marty and Tolstikhin, 
1998, Salters and Stracke, 2004], These source mantle CO 2 concentrations are inferred by 
four methods: (%) Cartigny et al. [2008] start with Nb concentration in MORBs, use Nb/C 
ratios, and assume an average degree of melting to calculate a CO 2 concentration in the source 
mantle, (ii) Saal et al. [2002] use MORB melt inclusions to measure CO 2 concentrations in the 
erupting mantle immediately prior to degassing, then assume an average degree of melting to 
calculate a CO 2 concentration in the source mantle. (Hi) Marty and Tolstikhin [1998] use 3 He 
concentrations in the ocean to infer a 3 He efflux from MORs. Then they apply a He/C ratio to 
estimate carbon efflux from MORs, which is matched to a CO 2 concentration in the erupting 
mantle using a (completely degassed) crustal formation rate of 21 km 3 /yr. CO 2 in the source 
mantle is then calculated by assuming an average degree of melting to generate MORB. (iv) 
Salters and Stracke [2004] start from major elements in the mantle and use a chain of element 
ratios to derive a CO 2 concentration in the depleted mantle. For all these approaches, CO 2 
concentration in the mantle is derived using several assumptions and has large uncertainties. 
The global MOR efflux of CO 2 is calculated using fewer assumptions and has less uncertainty. 
An alternative constraint on mantle CO 2 content is to match the global baseline emissions of 
the model to the estimated global MOR efflux by choosing a concentration of CO 2 in the source 
mantle of 215 ppmw (section 3.5). This is effectively substituting our average degree of melting 
into other’s calculations, and gives a result slightly higher than prior estimates of Cqo 2 - Fc> r 
simplicity, we use a generally accepted 125 ppmw Cqo 2 throughout this paper, but restate key 
results for the 215 ppmw value. This changes CO 2 emissions by a factor of 1.7, as emissions 
scale linearly with Cq,o 2 (see eqn.(3)). 

Subsequent sections present scenarios of changing sea-level. In these sections, plots depict 
emissions in terms of percentage difference from baseline emissions for the appropriate half¬ 
spreading rate; these percentage values apply to all highly incompatible elements and for any 
Cco 2 - Percentage results can be converted to CO 2 mass via figure 2. The deviations from the 
baseline emission rate are a consequence of non-zero z m in equation (3). By definition, total 
emissions are equal to this deviation summed with baseline emissions. 



Figure 2: (a) Ridge CO 2 emissions for constant sea level, and (b) focusing width Xf and 
full width x w of the melting region; all shown for varying half-spreading rate. Focusing 
width is equal to the width of the melting region when crustal thickness is less than 7 km 
(Uq < 1 cm/yr), and is otherwise limited such that crustal thickness does not exceed 7 km. 
The switch between these behaviours is marked by the grey dotted line at 1 cm/yr. Eqo 2 is 
computed with CO 2 concentration in the mantle of 125 ppmw. 

3.2 Single Sudden Change in Sea Level 

Imposition of a sudden sea-level change exposes the behaviour of the model. We use a steep, 
linear ramp, which gives a box function in S. Figure 3 shows the predicted MOR CO 2 emission 
rate Eqo 2 resulting from this sea level forcing. CO 2 emissions remain at the baseline level for 
~90 kyrs after the change in sea level and then there is a sharp rise in Eqo 2 , representing an 
8% increase in emission rate. This delayed response is due to the travel time of CO 2 from the 
base of the melting region to the ridge. After the 8% peak, E C q 2 falls sharply to about 3%, 
followed by a slow decay until 130 kyrs, after which there is a linear drop back to the baseline 
level over the duration of the box-pulse in S. The origin of this emissions pattern in figure 3(c) 
is explored in figure 4. 

Figure 4 demonstrates how the shape of the emission response curve in figure 3(c) is a 
consequence of the travel time r and its variation with respect to distance x from the ridge 
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Time (kyrs) 

Figure 3: Single step change in sea level. Plots show (a) sea-level change, (b) the negative 
rate of sea-level change, and (c) CO 2 emissions from a section of mid-ocean ridge, measured 
as percentage change from baseline emissions. The plate half-spreading rate and the per¬ 
meability constant at 1% porosity are Uq and K 0 , respectively. The three times marked on 
the plots correspond to (a—i,ii,iii) in figure 4. Negative S is plotted so that peaks in S and 
consequent peaks in Eqo 2 point in the same direction. 


axis. This travel time, shown in figure 4(c), ranges between r m ; n for melts originating in the 
mantle beneath the ridge axis and r max for distal melts that are focused laterally. The mixture 
of melts that arrives at the ridge at time t contains COo that was transported in melts that 
initiated along the base of the melting region at all x < Xf. These deep melts formed and 
began to segregate from their source at times in the past t — t(x). The CO 2 content of the 
segregated melt is different to the baseline case according to the S value at t — r{x). Therefore, 
we can calculate the deviation from the baseline by considering the rate of sea-level change 
(alternatively, z m ) acting on the melting region when each element of melt first segregated. 
This is represented in figure 4(c), and the integral of this plot with respect to x is directly 
proportional to E C q 2 in figure 3. 

The Eqo 2 response is clarified by again considering a sharp, linear drop of sea level that 
occurs over the interval 0 < t < 10 kyrs. The drop in sea level causes the depth of the base 
of the melting region to increase, importing additional CO 2 into the melting regime. In panels 
(a-i), (a-ii), and (a-iii) of fig. 4 we see the box pulse receding into the past as time progresses 
over t\ < t 2 < t 3 (figure 3c shows E C q 2 at these times). The start and end times of the pulse 
are projected onto t(x) in panel (b). At t\, the projection lines do not intersect r(x), meaning 
that the CO 2 perturbation generated by S has not yet reached the ridge axis. Therefore the 
emissions curve in figure 3(c) remains at the baseline level. At t 2 , the projection lines span 
r min ; the shallow slope of r{x) near r min means that C0 2 from a broad (30 km) region affected 
by the S pulse is arriving at the ridge axis at t 2 . This causes the spike in emissions shown in 
figure 3(c). At t 3 , the interval of sea-level change has receded far into the past. The only C0 2 
perturbation in melts arriving at the ridge is in distal melts from the base of the melting region 
at 50 to 58 km off-axis. The narrowness of this band translates to a reduced (but non-zero) 
value of Eqq 2 at t 3 in fig. 3(c). As the sea-level drop recedes further into the past, the emission 
rate drops to zero because the very distal melts (from x > xf) are not focused to the ridge axis. 
This process is animated in a video in the online supplementary materials. 
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Figure 4: The effect of a linear change in sea level on ridge CO 2 emissions. The three points 
in time shown in black, red, and green show the state of the system 65, 95, and 125 kyrs after 
the pulse began. Plots show: (a—i,ii,iii) rate of sea-level change from present into the past, 
for the three points in time marked in figure 3. (b) Travel time of melt from the base of the 
melting region to the ridge, increasing downwards. Travel time at x = 0 is slightly greater 
than 7 ~min because the sharp increase in lithospheric thickness within a few km of the ridge 
increases the column height; the effect of this increase in height exceeds the effect of faster 
on-axis upwelling. (c) Rate-of-change of the depth of first melting that acted on the melt 
currently arriving at the ridge axis, when that melt began to segregate. Dashed lines mark 
how the travel time converts S(t) to z m (x). The integral of a coloured line from 0 to Xf in 
panel (c) is directly proportional to Eqo 2 a t the corresponding time in figure 3(c). This is 
expressed by the integral of z m with respect to x in equation (3). A video animating this plot 
and Eq o 2 ° ver time is included in the online supplementary materials. 


In the limit of vanishing duration of sea-level change, the pulse of S approximates to a 
Dirac delta function. Hence the Eqo 2 response shown in figure 3(c) is an approximation of the 
Green’s function. Conceptually, the emissions response for any sea-level time series could be 
approximated by convolution with a Green’s function approximation like figure 3(c). 

From this discussion, it is clear that changes to t(x) will alter how MOR emissions respond 
to changing sea-level. Equation (7) for travel time shows that plate half-spreading rate Uq and 
permeability constant Kq are the key parameters affecting the melt travel-time r. Higher values 
of either lead to shorter travel times, giving a higher peak in Eqq 2 over a shorter time period, 
with this peak occurring sooner after the causative sea-level change. This is explored in more 
detail in the following section. 

3.3 Oscillating Sea Level 

We now consider oscillatory sea level and discuss the concepts of lag and admittance. 

Figure 5 shows a pair of oscillating sea level scenarios and their predicted Eqq 2 variation. 
The left column (i) shows a time series of alternating box-pulses in S] the right column (ii) shows 
a sinusoidal sea-level variation. Figure 5(c-i) has an oscillating series of peaks in Eqo 2 resulting 
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Figure 5: Sawtooth and sine waves in sea level starting at t = 0. Plots show (a) sea-level, 
(b) negative rate of sea-level change, (c) CO 2 emissions from a section of mid-ocean ridge. 
Considering the left column, each box-pulse in S produces an emissions peak/trough as in 
figure 3(c) with interference where they overlap, giving the net result seen in (c-i). The steady 
upward shift in the -Eco 2 -P ea ks until ~150 kyr is due to this overlap. If the first box-pulse in 
S had been in the opposite direction, Eco 2 peaks in (c) would be mirrored across Eqo 2 — 0. 


from a series of S box-pulses in figure 5(b-i). This S series is equivalent to summing single 
box-pulses from figure 3(b), with suitable offset and amplitude. Similarly, the CO 2 emissions 
can be represented as a sum of offset emissions spikes from single, linear changes in sea level. 
The train of emissions peaks in figure 5 (c-i) and (c-ii) stabilises in amplitude after t « r max ; this 
transient represents the spin-up time of the model, associated with the tail of excess emissions 
shown in fig. 3. 



Period (kyr) U 0 (cm/yr) log(/\ 0 ) (m 2 ) 

Figure 6: Lag and travel time for varying (a) sinusoidal sea-level period, (b) half-spreading 
rate, and (c) permeability constant. The lag is calculated as the time between a peak in —S 
and the corresponding peak in Eqo 2 f° r sinusoidal sea level. 

Figure 5(ii) shows sinusoidal sea level and provides the context to define the lag metric. Lag 
£, is the time between a peak in S and the corresponding peak in emissions. Because the time 
interval around r min kyrs before t has the largest influence on Eco 2 at t, the Eco 2 signal should 
lag S by about r min . However, the lag will not be exactly r min , as the influence of S on E C q 2 is 
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felt up to r max years after the change in sea level. Thus the exact value of lag will be slightly 
greater than r m j n and we expect this difference to depend on the period of sea-level oscillation 
relative to r max —r m | n . In particular, when the period approaches and exceeds 2(r max — r m i n ), the 
lag becomes equal to the mean melt-transport time r mean . Figure 6 shows lag, r min , T mean , and 
r max for varying half-spreading rate, permeability constant, and sinusoidal sea-level period. We 
note that r min < C < r mean <C r max and therefore we assume C ~ r min with a small, systematic 
error. 


Uo fixed = 3 cm/yr 



A'o fixed = 10 12 m 2 




Sea level period (kyr) Sea level period (kyr) 


Figure 7: Absolute and relative admittance for varying half-spreading rate, permeability, 
and sinusoidal sea-level period. Plots hold either Uo or K 0 fixed. The magnitude of sea-level 
change is constant for all periods. 


Sinusoidal variation of sea level also provides a context in which compute admittance. Ad¬ 
mittance is the ratio of the response amplitude to the forcing amplitude as a function of the 
sinusoidal forcing period. We define two versions of admittance: absolute admittance, with 
units of kilograms CO 2 per metre of ridge per year per 100 m of sea-level change, and relative 
admittance, with units of percentage change from baseline emission rate per 100 m of sea-level 
change. The latter is the absolute admittance divided by the baseline emission rate. Figure 7 
shows absolute and relative admittance and how they vary with changing sea-level period, half¬ 
spreading rate, and permeability constant. We discuss both the trends and the oscillations of 
these curves, starting with absolute admittance. 

Absolute admittance (panels (a) and (b) of figure 7) depends on the period of sea-level 
oscillation, the permeability, and to a lesser extent, on the half-spreading rate. We consider 
these in turn. Shorter periods of sea-level variation at constant amplitude give larger values 
of (S’! and |i m |, and hence increase the temporal variation of /co 2 - This causes increased 
deviation of ridge emissions from baseline. Increased mantle permeability and spreading rate 
both reduce the melt travel-time from the base of the melting region. In the melt travel-time 
model of equation (7), the permeability appears as K 0 , while the spreading rate Uo is directly 
proportional to the mantle upwelling rate W m . A reduced melt travel-time implies a smaller 
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difference between r max and r m i n , and therefore a focusing of CO 2 from the base of the melting 
region to the ridge axis over a shorter interval in time. The uncertainty in mantle permeability 
translates to a broader spread of the absolute admittance curves than does the range of half¬ 
spreading rates considered. 

Relative admittance is equal to the absolute admittance normalised by the baseline emissions 
rate. The baseline depends on half-spreading rate but not on permeability (fig. 2(a)). We 
therefore see a difference between absolute and relative admittance in figure 7(b) and (d). For 
slow-spreading ridges, which have a low baseline emissions rate, the normalised variance (and 
hence the relative admittance) is larger. 

The oscillations superimposed on the primary admittance trend are not physically signifi¬ 
cant, but are readily explained. They arise from variation in the number of sea-level half-cycles 
that fit into the time interval from t — r max to t — r m this is the time interval over which 
S(t — r) can contribute CO 2 emissions at time t. For oscillatory sea level, each positive or 
negative peak in S has an opposing emissions effect relative to the prior negative or positive 
peak. If there is an unmatched peak affecting the bottom of the melting region, the amplitude 
of Eqq 2 variations is larger and the admittance is higher. 

Broadly, the patterns of admittance and lag imply that, in terms of CO 2 emission variation, 
the dominant sea-level changes will be those with large amplitude and short period changes. 
The emissions variation associated with such changes will lag the forcing by approximately T m \ n . 
The modelled magnitude and lag of Eqq 2 changes are affected by both K 0 and U 0 . 

3.4 Reconstructed Pleistocene Sea Level 

The simple scenarios of sea-level variation presented above give insight into the behaviour of 
the model, but are not representative of the variations that have occurred naturally, over the 
past million years. We move, therefore, to a model forced by the time-series of reconstructed 
global sea level from Siddall et al. [2010], shown in figure 8(a). Other reconstructions exist, but 
the differences between them are small enough that we follow Crowley et al. [2015] and consider 
only this one. Siddall et al. [2010] record data every 3 kyrs and, based on their reconstruction, 
the highest rates of sea-level change (fig. 8(b)) meet the condition rnax(i m (t)) < W m required 
for validity of equation (3). 

Figure 8(c) shows the result of applying reconstructed sea level to ridge E C q 2 . There is a 
±8 % range in Eqo 2 for moderate half-spreading rate and permeability. The Eqo 2 curve is, 
qualitatively, an offset version of S with small variations smoothed out — as expected from 
Eqo 2 being approximated by convolving S with the emissions response in figure 3(c). Within 
this framework, we now consider how to apply the model to global MOR emissions. 

3.5 Global Mid-Ocean Ridges 

The global MOR system is composed of ridge segments spreading at different rates, ranging 
from the ultra-slow Gakkel ridge to the fast East Pacific Rise. The baseline emissions depend 
on the half-spreading rate, as does the character of the emissions response to sea-level change. 
The global response to sea-level change should therefore be computed as the segment-scale 
response, integrated over the global MOR system, 

/"^MOR 

Gco 2 (t,Ko) = / E C o 2 (t, K 0 , U 0 (l))dl, (9) 

where Gco 2 is global MOR emissions of C0 2 in kg/yr, l is arc length along the ridge, and Tmor 
is the total length of the MOR system. This integral can be approximated by discretising the 
half-spreading rate into bins [/ 0 ? and summing the local response in each bin. A weighting is 
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Time (kyrs) 


Figure 8: Reconstructed sea level. Plots show (a) sea-level change, (b) the negative rate 
of sea-level change, and (c) CO 2 emissions from a section of mid-ocean ridge. The lag is 
90 kyrs. 


applied to each entry in the sum to account for the total length of segments with half-spreading 
rates in that bin. The sum is written as 

N 

Gco 2 (t,K 0 ) = Y E C q 2 ( t , K 0 , U 0i ) Li(U 0 i), (10) 

i= 1 

where N is the total number of spreading rate intervals to sum over and L,(t/oi) is the total 
length of MOR in a particular spreading-rate bin. The local emission rate Eqo 2 in each bin is 
computed by adopting the average half-spreading rate of the bin and assuming that sea-level 
change is eustatic — the same for all segments globally. 

Gale et al. [2013] provide a catalogue of segment lengths and spreading rates for the global 
MOR system; the total ridge length is 61000 km with a mean half-spreading rate of 2.5 cm/yr. 
A histogram of these data is plotted in figure 9(a). Figure 9(b) shows, for each spreading- 
rate bin, the rate of baseline emissions per metre of ridge. Panel (c) then shows the product 
of ridge length and emissions rate per metre, giving the total emissions rate associated with 
each spreading rate bin. These are summed in accordance with eqn. (10) to give the total 
baseline global response. The global baseline emission rate thus predicted is 53 Mt C0 2 per 
year assuming a sub-ridge mantle C0 2 concentration of 125 ppmw. This can be compared to 
other estimates of 91 ± 45 MtC0 2 /yr [Coltice et ah, 2004, Cartigny et ah, 2008, Marty and 
Tolstikhin, 1998]. If we instead constrain the model to have baseline emissions of 91 MtC0 2 /yr, 
it requires a sub-ridge mantle C0 2 concentration of 215 ppmw (0.65 kg per m 3 ). 

Before applying the reconstructed sea-level forcing to the weighted global emissions sum 
in eqn. (10), we consider the simple sea level forcing that was used to probe the behaviour 
of E C o 2 hr fig- 3. This linear ramp in sea level is applied to compute the global emissions 
response in figure 10(a) for a range of mantle permeabilities. Global emissions in figure 10(a) 
are, unsurprisingly, more complex than the Eqo 2 equivalent, as they consist of a summation 
of N Eqq 2 peaks from figure 3, weighted according to the ridge length in each bin and offset 
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Figure 9: Quantities associated with discrete bins of half-spreading rate, (a) The total 
length of MOR [Gale et ah, 2013]. (b) Baseline emissions per metre of MOR from our 
model, (c) Baseline global CO 2 emissions as the product of quantities in panels (a) and (b). 
Calculated for 0.38 kg CO 2 per m 3 of upwelling mantle (125 ppmw). Baseline emissions are 
independent of K 0 . 


due to the variation in r with f7 0 n Compared with figure 3, Gco 2 l ias a smaller percentage 
difference in the rate of CO 2 emissions and is spread out over a longer time-period. 

Finally, we apply the reconstructed sea-level time series to the global model. Figure 10 shows 
Gco 2 f° r reconstructed sea level, calculated and plotted for a set of three values of permeability 
constant K 0 in panels (d-i), (d-ii), and (d-iii). These curves demonstrate a reduction in Gco 2 
range from 6.6 MtC02/yr for K 0 = 10 -11 m 2 to 1.2 MtC02/yr for K 0 = 10 -13 m 2 . The 
reduction in Gco 2 range occurs because K 0 affects the range of r min that is implicit in the 
global sum in eqn. (10). A large value of K 0 gives higher overall permeability, shorter melt- 
transport times, and a global range in r min that is smaller. Therefore the emissions response to 
a box-pulse in S (fig. 10(a)) is temporally concentrated and attains a higher peak value. Hence, 
larger K 0 causes greater amplitude of variation in Gco 2 for reconstructed sea level. 


16 


































-900 -800 -700 -600 -500 -400 -300 -200 -100 

Time (kyrs) 


Figure 10: Global MOR emissions. Panel (a) shows CO 2 emissions from the global MOR 
system for a step change in sea-level identical to that of fig. 3(a). Subsequent panels consider 
reconstructed sea level and its effect, (b) Reconstructed sea level; (c) the negative rate of 
sea-level change; and (d-i—iii) CO 2 emissions from the global MOR system. The sea level 
time-series has been applied further back in time than shown, such that the left-most point 
in f?co 2 is affected by more than r max kyrs of prior sea-level change. For (d-i,ii,iii) the lags in 
Gco 2 are, respectively, ~60 kyrs, ~120 kyrs, and ~250 kyrs. 

4 Discussion 

The model prediction for global MOR CO 2 emissions at constant sea level using a mantle CO 2 
concentration of 125 ppmw is 53 Mt CO 2 per year (14 Mt carbon per year), a value that 
is within error of estimates from analyses of mid-ocean ridge basalts of 91 ± 45 MtC02/yr 
[Coltice et ah, 2004, Cartigny et al., 2008, Marty and Tolstikhin, 1998]. However, note that our 
prediction is not independent of MORB studies because our choice of Cco 2 is based on observed 
MORB chemistry. If we instead choose to constrain the model to fit published CO 2 emissions, 
we require a source mantle CO 2 concentration of ~215 ppmw. 

The ranges in global CO 2 emissions under reconstructed variation of past sea level (fig. 10 
for 125 ppmw CO 2 ) are 1 MtCC^/yr, 3 MtCC^/yr, and 7 MtCC^/yr for mantle permeabilities 
at 1% porosity of 10 -13 m 2 , 10~ 12 m 2 , and 10~ n m 2 respectively. These idealised predictions 
assume that 100% of C0 2 transported to the ridge axis is degassed into the oceans, rather 
than retained in the crust or mantle. This may, in fact, be rather accurate; Cartigny et ah 
[2008] estimated that over 80% of CO 2 in primitive MORB is degassed near the ridge axis. 
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Furthermore, 100% degassing is assumed in the papers that calculate the source mantle CO 2 
concentration we use. It has previously been speculated that variations in ridge CO 2 emissions 
would have an effect on global climate [Huybers and Langmuir, 2009, Lund and Asirnow, 2011]; 
it is beyond the scope of the present manuscript to investigate this question, though it is a 
target for future work. 

Uncertainty in mantle permeability translates to uncertainty in both the amplitude of vari¬ 
ations and the lag of global CO 2 emissions from MORs. There are various experimental con¬ 
straints on mantle permeability [e.g. Miller et ah, 2014, Connolly et al., 2009]; these tend to 
agree on the scaling with porosity, but disagree on the magnitude of K 0 . Furthermore, per¬ 
meability is sensitive to grain size, a parameter that is poorly known for the mantle beneath 
MORs [although see Turner et ah, 2015]. Our chosen range of K 0 is intended to accommo¬ 
date these uncertainties, as well as represent an effective permeability for melt transport that 
may be channelised into high-porosity, high-permeability dunite channels [e.g. Jull et ah, 2002, 
Kelemen et ah, 1995]. Our Kq range of 10~ 13 -10~ n m 2 encompasses a change in the amplitude 
of Gco 2 variation by a factor of 5, a difference in lag of 200 kyrs, and qualitative difference 
in the time-series of Gco 2 (fig- 10)- Therefore K 0 represents a leading source of uncertainty 
in the model. Uranium series disequilibria may provide an independent constraint on magma 
travel time from the base of the melting region [e.g., Jull et ah, 2002], although interpreting the 
various species in the decay chain is fraught with complexity. Preservation of 230 Th disequilib¬ 
rium (half-life of 75 kyrs) suggests a permeability of K 0 > 1CD 12 m 2 , and community consensus 
similarly favours K 0 at the higher end of our considered range. 

Our model is based on the assumption that melt travels vertically from the depth of first 
melting to the top of the melting region and is then focused along a sloping, high porosity 
decompaction channel to the ridge axis [Sparks and Parmentier, 1991]. The travel time of the 
vertical flow is modelled by a ID compaction column; we assume that transport in the decom¬ 
paction channel is instantaneous. The systematic error introduced by the latter assumption is 
zero on the ridge axis (x = 0), and increases with distance x from the ridge axis. This means r 
plotted in figure 4(b) is more accurate at small x , but increasingly underestimates r for larger 
x. Therefore, assuming r min is accurately modelled, r max is probably too small, such that Eq o 2 
in figure 3 should have a longer tail on the right of the graph. However, long tails have little 
effect on the resultant Eqo 2 pattern for complex sea-level changes or on the Gco 2 pattern for 
reconstructed sea level. Therefore the overall effect of including a finite travel time along the 
high porosity channel would be to make a small adjustment to the Eqo 2 response. This suggests 
that assuming instantaneous travel time along the channel has little effect on the results of the 
model. 

Another assumption made is that travel time is constant with respect to time, despite 
changes in melting rate (and thus, porosity) caused by changing sea-level. This follows the 
approach in the reduced model of Crowley et al. [2015], where the perturbations in porosity 
were taken as negligible disturbances to the travel time in a steady-state compaction column 
solution from Hewitt [2010]. 

A more significant assumption underpinning the model is that carbon dioxide behaves as a 
perfectly incompatible tracer — meaning that carbon concentration does not affect the mantle’s 
physical or thermodynamic properties. However, carbon is not a trace element. In contrast, 
experiments by Kono et al. [2014] document the very low viscosity of insipient, carbon-rich melts 
present at small melt fraction below the base of the silicate melting region. The experiments 
also show that viscosity rises sharply as the carbon is diluted by silicate melting. It would be 
challenging to capture this variability in models, especially since the wetting properties (and 
hence the mobility) of carbon-rich melts are poorly constrained. 

A more significant concern, however, is that treating volatiles as trace elements neglects 
their thermodynamic effect on melting. Small mantle concentrations of carbon affect the depth 


18 


at which melting begins, though the melt fractions produced by this incipient melting are 
probably less than a few tenths of a percent [Dasgupta et ah, 2013]. Our model assumes that 
these melts do not segregate until the onset of silicate melting. At such small porosity, it is 
unclear whether these carbonated melts can percolate. However, water-induced melting at the 
wet peridotite solidus of ~90-120 km [Asimow and Langmuir, 2003, Dasgupta et ah, 2013] 
increases the melt fraction. Again, the threshold of interconnectivity for such melts is not 
known, so it is possible that such deep, hydrous melts do not segregate, or do so very slowly. If 
the 230 Th disequilibrium observed in young MOR lavas originates with melt segregation in the 
presence of garnet [Stracke et al., 2006], it would support the hypothesis of efficient segregation 
of hydrous melts (although other hypotheses also fit the observations). Overall, our model 
depends only on the presence of a pressure-dependent boundary that separates non-segregating 
melts, below, from segregating melts, above. The sharpness that is required of this boundary 
is unclear. 

Finally, our model assumes a chemically and thermally homogenous mantle, which is cer¬ 
tainly not true of the natural system [ e.g ., Dalton et ah, 2014], No data exists that would allow 
us to accurately incorporate small-scale (< 100 km) heterogeneity in model. If such heterogene¬ 
ity is pervasive and at scales of ~10 km or smaller, it would affect the style of melt transport 
[Katz and Weatherley, 2012], with fertile regions creating pathways for rapid melt transport 
through the melting region. It may be the case that this is captured by a high effective per¬ 
meability, though this is probably not a testable hypothesis. Large-scale heterogeneity would 
leave the mantle homogenous over the scale at which we calculate Eqq 2 , so the underlying melt 
transport model would be unaffected, though parameters would need to be adjusted according 
to the oceanic region. It seems likely that such variations would cancel in the integral for global 
CO 2 emissions rate. 

Previous authors have suggested that changing sea level might affect MOR CO 2 emissions 
almost instantaneously [Huybers and Langmuir, 2009, Lund and Asimow, 2011, Tolstoy, 2015]. 
Their assumption is that pressure changes instantaneously affect melting rates, MOR volcanic 
productivity and therefore, assuming constant CO 2 concentration in the erupted melt, MOR 
CO 2 emissions. We disagree with the assumption that CO 2 concentrations would be constant. 
CO 2 is highly incompatible and therefore additional melting acts to dilute the constant mass 
of CO 2 in the melt. However, after including these effects, our model can calculate whether 
changing sea level affects MOR CO 2 emissions (see appendix A. 4). To leading order there is no 
effect; the reduced (or increased) concentration of CO 2 in the melt counteracts the increased 
(or reduced) rate of melt delivery to the ridge axis. 

We would like to be able to compare the model to data, but there is no dataset of global 
MOR CO 2 flux over time. Atmospheric CO 2 concentration from Antarctic ice cores [Bereiter 
et ah, 2015] is an existing dataset that might record some influence of Gco 2 - However, there 
are many strong, nonlinear controls on atmospheric CO 2 and the relationship between Gco 2 
and atmospheric CO 2 will not necessarily be constant over time. Nonetheless, comparing 
atmospheric CO 2 to Gco 2 does offer at least one useful insight: if there is no correlation between 
these quantities within the model’s reasonable parameter space, it would be an indication that 
MOR CO 2 emissions have no effect on atmospheric CO 2 . We find a correlation in appendix A.5. 
The best correlation, without considering the very significant complicating factors associated 
with climatic feedbacks, occurs at K a = 10~ 1L37 with a lag of 73 kyrs. This correlation does 
not prove any causative link between Gco 2 and atmospheric CO 2 . 

Variable MOR CO 2 emissions’ effect on atmospheric CO 2 will vary over time, as the fraction 
of MOR CO 2 emissions into the intermediate ocean that reach the atmosphere is not constant. 
This fraction depends upon an array of factors that affect ocean alkalinity, ocean upwelling 
patterns and CO 2 concentration in the atmosphere, and these factors vary on sub-glacial/glacial 
timescales (as do numerous processes interacting with them, e.g., the biosphere) [Broecker, 1982, 
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Walker, 1986, Brovkin et al., 2012, Yu et al., 2014]. This would have to be taken into account 
if inferring any climate implications from our model. 

The MOR carbon flux may vary over time in other ways that we have not considered here. 
For instance, it is plausible that the intensity of hydrothermal circulation varies with sea level, 
driven by variations in melt supply [Lund and Asirnow, 2011, Crowley et ah, 2015]. If this 
is the case, hydrothermal variations would have a different lag than that of CO 2 emissions. 
Hydrothermal systems have been proposed as both a CO 2 sink, with hot seawater transforming 
basalts to clay [Gysi and Stefansson, 2012], and a CO 2 source, with hydrothermal fluids trans¬ 
porting CO 2 from magma to the ocean [Sakai et ah, 1990]. The rate of both these processes 
might scale with hydrothermal circulation, although it is not clear whether the net effect would 
be to increase or decrease MOR CO 2 emissions. 

5 Summary 

The model presented above builds on the reduced model of Crowley et al. [2015] to calculate 
the efflux of a highly incompatible chemical component from a mid-ocean ridge, and how that 
efflux would vary with changes in sea level. It is based on a description of melt transport 
through a homogenous mantle and assumes perfect incompatibility of the component. This 
leads to a simple but physically consistent model of chemical transport through the melting 
region beneath a ridge. The model calculates total melt supply rate and global background 
emissions of CO 2 that are consistent with data and prior estimates. 

In the model, changing sea level affects the depth of first silicate melting; this alters the 
rate at which CO 2 enters the melting region, segregates from its mantle source, and (some time 
later) arrives at the ridge axis and is degassed into the ocean. 

The MOR emission rate of C0 2 is predicted to vary by up to 12.5% when the model 
is forced with reconstructed Pleistocene sea level variation (7 MtC0 2 /yr for 125 ppmw Cqo 2 ; 

11 MtC0 2 /yr for 215 ppmw Cqo 2 )• There is uncertainty in the predicted magnitude and timing 
(relative to sea-level forcing) of this effect, as two parameters of the model — Cco 2 and Ko 
- are weakly constrained by existing data. However, within reasonable ranges of the model 
parameters, the amplitude of global MOR CO 2 emission-rate variation will remain on the order 
of several MtC02/yr. The total difference in the mass of CO 2 emitted from MORs during 
sea-level driven deviations from global baseline CO 2 is up to ~80 Gt CO 2 for high permeability 
and 125 ppmw Cco 2 ( see figure 10(d-i) ). This is 4% of the pre-industrial CO 2 mass in the 
atmosphere of 2190 Gt, or 0.06% of pre-industrial CO 2 in the oceans [IPCC: Solomon et ah, 
2007, fig 7.3], 

Our results indicate that the CO 2 emissions from mid-ocean ridges are temporally variable in 
response to sea-level change. These results align with the hypothesis of Huybers and Langmuir 
[2009]; however, whereas they assumed an immediate emissions response, we show that the mid¬ 
ocean ridge CO 2 emissions response lags the sea-level forcing by approximately the minimum 
travel time of CO 2 through the melting region — at least 60 kyrs. Therefore MOR CO 2 
emissions cannot feed back into the sea-level change that caused them; instead these CO 2 
emissions will enter the climate system during the next glacial cycle. Further work is needed 
to assess the climatic impact of the variable CO 2 emissions predicted for mid-ocean ridges in 
this paper and for subaerial volcanoes by Huybers and Langmuir [2009]. 
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A Appendix 

A.l Depth of First Silicate Melting 

The depth of first silicate melting occurs at the intersection of the local mantle temperature, 
which is adiabatic prior to first melting, and solidus temperature. We take a standard adiabatic 
temperature expression and approximate it with a Taylor series truncated at first order in 
pressure to obtain 

T { Z ) = f-^(pz-p w S) , (11) 

pc 

where symbols are as defined in the main text and we have added a pressure term for varying 
sea level. The mass of a column of water is p w S, provided the thermal expansion of the ocean is 
much less than S. Fortunately, oceanic thermal expansion causes less than 1% of total glacial 
sea-level change [McKay et ah, 2011]. Freshwater density is used for p w , as salt remains in 
the ocean when water is removed, and ocean water inputs are fresh. Coordinates are upwards- 
positive with origin at the MOR, therefore z is negative and increasing sea level S is positive. 
We model the solidus temperature following Katz [2008] and Hewitt [2010], using a solidus that 
is linear in pressure and a single composition parameter. We consider pressure terms due to 
the mass of mantle above the solidus, and sea level deviations from reference conditions. This 
gives 

Solidus = T Sie{ - jg (.pz - p w S ) . (12) 

We ignore bulk ocean+atmospheric weight as this is accounted for in Ts ief . Setting Tsoiidus in 
eqn. (12) equal to T(z) in eqn. (11) and solving for depth we obtain eqn. (5) in the main text. 


A.2 Upper Boundary of the Melting Region 

The upper boundary of the melting region, like the lower boundary, is a region where the solidus 
temperature matches the local mantle temperature. Local mantle temperature in the melting 
region can be modelled as the adiabatic temperature profile (eqn. (11)) plus terms for (i) ther¬ 
mal energy lost to latent heat during melting of (z — z m )IIL/pc, and (ii) conduction near the 
cold ceiling (i.e., the ocean floor) from the half-space cooling model. This approach assumes 
superposition of temperature fields without cross terms correcting for deviations from the as¬ 
sumptions of each model (e.g., half-space cooling assumes a constant background temperature 
with respect to depth). Matching this to solidus temperature from equation (12) gives 


Fsvef - 7 pgz 



n l 

pc 


(z - z m ) - (T - T 0 ) erfc 




(13) 


where L is the latent heat capacity of the mantle. 

Equation (13) cannot be rearranged to give z in terms of x. However we can get x in terms 
of z, as shown in equation (14), an expression shown to accurately match melting region from 
numerical studies in Crowley et al. [2015]. This justifies the superposition of temperature fields 
assumption and gives us a reference to benchmark against, 


xi(z) 


Upz 2 

4k 



7 pg - agT/c-UL/pc 
T-Tn 


(z - Z m ) 


-2 


(14) 


To get z in terms of x, we discard all depth-dependent temperature terms in equation (13) 
except half-space cooling. This approach is equivalent to plotting an isotherm in a half space 
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cooling model. Thus 


erfc " 1 (?T^) ' (15) 

Figure 11 compares equation (15) to the benchmark (eqn. (14)). For a lithospheric boundary 
temperature of ~1550 K — close to the solidus temperature of mantle — equation (15) is a poor 
approximation of the real melting region. However, unphysically high boundary temperatures 
of > 1640 K give good approximations of the melting region, as this compensates for the 
temperature effects not explicitly accounted for. Therefore, we use equation (15) with T) = 
1643 K to approximate the upper boundary of the melting region. 


a 

a 
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Figure 11: Melting-region calculation from [Crowley et ah, 2015], compared to isotherms 
of the half-space cooling solution for 1565 K (a realistic solidus temperature) and 1643 K 
(a best fit isotherm). The 1643 K isotherm is used as the upper boundary of the melting 
region for this paper. All curves are calculated for a mantle potential temperature of 1648 K, 
surface temperature of 273 K and a half-spreading rate of 3 cm/yr. Plotted quantities all 
scale similarly with Uq and thus isotherms have the same relative error across all spreading 
rates. 


A.3 Width of the Melting Region 

Equation (14) shows a finite width of the melting region defined by dxi/dz = 0 (see figure 11 
at £ ~ 55 km, x ~ 130 km). The derivative of equation (14) is 


d xi(z) 
d2 


U 0 z (2 erfc -1 {B(z - z m )) - 


Ak 


where: B = 


[erfc l {B{z - Zm))] 

IP 9 ~ agT/c-UL/pc 


T-Tn 


(16) 

(17) 


We have not found an analytical solution to equation (16) for dxi/dz = 0. However, 
computational investigation of equation (14) for varying Uq shows that the depth of the widest 
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point of the melting region is constant. Calculating xi at this fixed depth gives the maximum 
width of the melting region. Thus, 


max(xj) = x w = (erfc 1 [B (. z w - z m )}) 


-2 


(18) 


where z w is the depth of the widest point of the melting region. 


A.4 MOR C02 Emissions Response to Changing Melting Rate 

We wish to know the effect of changing melting rates (due to pressure changes) on the MOR 
efflux of CO 2 . Following the model used in the main text, the CO 2 emissions rate from a 
one-dimensional column in the melting region is: 


Bco 2 ~ 


RC 


melt 

CO2 


(19) 


where R is the rate at which melt upwells per unit cross-sectional area (kg m -2 s” 1 ), and Cqq* 
is the mass concentration of CO 2 in the melt, both defined at the top of the column. The right 
hand side of eqn (19) is affected by changing pressure and its consequent effects on melting 
rate. 

Increased melting rate has two effects. Firstly, it increases porosity, which induces faster 
melt upwelling (increasing R). Secondly, it dilutes C0 2 in the melt, reducing C™g 1 * (C0 2 is 
highly incompatible, thus its mass in the melt stays constant regardless of changes to the total 
melt mass). These two effects act in opposition, and the net result on Eco 2 is not immediately 
obvious. However, it can be evaluated by calculating the rate of change of Eqo 2 in response to 
changing pressure, 


d£co a dfi^ melt , dCS5>* 
df ■ di Lc °“ + di 


( 20 ) 


This calculation requires expressions for Cqq*, R and their derivatives. Taking CO 2 as perfectly 
incompatible, 


‘-'CO2 — 
j fmelt 

QU CQ 2 _ 

df 


£co 2 
F ’ 
£co 2 d-F 
F 2 df 


(21a) 

(21b) 


where F is the degree of melting. For a linearised solidus, following Katz [2008] Hewitt [2010], 
this is 


F=^-(P- P m ) 
P 9 

d F _ n p w dS 

df p 2 g df 


(22a) 

(22b) 


where P is pressure, P m is the pressure at first silicate melting and we have used P = p w S. 

R is the product of the local porosity 0 and melt upwelling velocity w. For a compaction 
column, following Hewitt [2010], these are 


f KpApg 

V Vf 

( KpApg 

V Vf 


n w„ 

p 2 

HR» 

p 2 


l-i 


(P ~ Pm) 


(P-Prr 


1 — 


(23a) 

(23b) 
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Thus R = 4>w is 


d R 
d t 


R 



(24b) 


(24a) 


Substituting equations (21a), (21b), (22b), (24a), and (24b) into equation (20) then rear¬ 
ranging gives 



(25) 


where, following eqn (22a), the term inside parentheses is 1 — F/F = 0. Therefore, to leading 
order, the instantaneous effect of changing sea-level on MOR CO 2 emissions (through changes 
to melting rate and eruption volume) is zero. Thus, suggestions that MOR CO 2 emissions 
should scale with eruptive volume are likely incorrect; it is only more compatible elements that 
increase emissions rates with increased melting. 

However, the argument presented above is a simplification. For instance it assumes all melt 
that reaches the ridge is erupted instantaneously. A more detailed model might consider magma 
reservoirs at the MOR which could be vented by abrupt changes in pressure. However, such 
ventings would only be fluctuations against the background supply of C0 2 to the ridge, which 
(as shown above) does not vary with melting rate. 

A.5 Correlation to Atmospheric C02 

Whilst there is no dataset of global MOR CO 2 flux over time, we can compare the model’s 
calculated MOR CO 2 flux to the related quantity of atmospheric CO 2 . Specifically, we can use 
atmospheric CO 2 concentrations from amalgamated Antarctic ice core data [Bereiter et ah, 
2015]. We do not believe there will be a simple correlation between global MOR carbon 
emissions Gco 2 an d atmospheric CO 2 , nor do we think that this carbon flux is large enough to 
dominate the atmospheric CO 2 variations. However, if there is no correlation between these 
quantities within the model’s reasonable parameter space, it would be a strong indication that 
MOR CO 2 emissions have no effect on atmospheric CO 2 . 

In plotting atmospheric CO 2 (mass) against Gco 2 (mass/yr), figure 12 assumes that the at¬ 
mosphere instantaneously adjusts its equilibrium CO 2 concentration in response to any change 
in outgassing at MORs. This is a reasonable timescale approximation, as the equilibrium 
timescales are short compared to our study. However, assuming a linear change in equilibrium 
atmospheric CO 2 concentration in response to Gco 2 , despite atmospheric CO 2 having multi¬ 
ple feedbacks and the time-varying fraction of MOR emissions entering the atmosphere (see 
section 4), is a considerable flaw. 

Figure 12 shows that a correlation between these series is possible within the reasonable 
parameter space of our model’s least-well constrained parameter K 0 , however this is not a proof 
of causal links. The process we have used could even be considered circular. Atmospheric CO 2 
and sea level are strongly correlated and our model is driven by (rate of change of) sea level; it 
is perhaps unsurprising that the model output is also correlated with atmospheric CO 2 . 
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Figure 12: Best fit of global MOR CO 2 emissions (125 ppmw Cco 2 ) to atmospheric CO 2 
mass. Atmospheric C0 2 mass was calculated using 1 ppmv C0 2 = 7.81 Gt C0 2 . Atmospheric 
C0 2 mass has been plotted to demonstrate that Gco 2 variations of a few Mt/yr (10 9 kg/yr) 
for 10’s of kyrs (10 4 yr) have less mass than the change in atmospheric COo mass (10 14 kg) by 
1-2 orders of magnitude. Therefore correlation between these quantities would have to be the 
result of feedbacks on atmospheric CO 2 . Best fit was found by varying Kq from ICG 13 to 10~ n 
and finding the largest Pearson product-moment correlation coefficient (0.45). Maximising 
Spearman (0.48) or Kendall (0.32) correlations produces the same Kq value. The lag for best 
fit K 0 of 10" 11 ' 37 m 2 is 73 kyrs. 
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